module data
   !
   !  SUBROUTINES
   !
   !  1) get_parameters
   !  2) write_parameters
   !  3) allocate_variables
   !  4) initialize_variables
   !  5) date_time
   !
   ! Last update: 19 Oct 2012
   !
   use portlib

   implicit none
   !
   ! Global variables
   !
   integer :: nx     ! Number of real + ficititious volumes in the csi direction
   integer :: ny     ! Number of real + ficititious volumes in the eta direction
   integer :: nxy    ! Total number of real + fictitious volumes
   integer :: kg     ! Kind of grid
   integer :: ctd    ! Kind of grid in axial direction
   integer :: coord  ! Coordinate system (1=cylindrical; 0=cartesian)
   integer :: specie ! chemical specie (1=H20; 2=AIR)
   integer :: modvis ! Viscosity model (0=Euler; 1=Navier-Stokes)
   integer :: ccTw   ! boundary condition (0=adiabatic; 1=prescribed temperature)
   integer :: ig     ! Number of volume (x direction) of the throttle
   integer :: ito    ! First iteraction (ito = 1 at first iteraction and other value for reload)
   integer :: itmax  ! Maximum number of iteractions of the time evolution
   integer :: ittol  ! number of iterations below the tolerance for stopping criterion
   integer :: imax   ! Maximum number of iteractions for the pressure correction
   integer :: nitm_u ! Maximum number of iteractions for the MSI method for u, v and T
   integer :: nitm_p ! Maximum number of iteractions for the MSI method for p
   integer :: reload ! Upload backup data and continue computation (0=no; 1=yes)
   integer :: wbkp   ! Frequency of saving backup data
   integer :: wlf    ! Frequency of printing in the listing file
   integer :: sem_a  ! open result files (0=yes; 1=no)
   integer :: sem_g  ! Visualize the plot (0=yes; 1=no)
   integer :: w_g    ! Frequency of writing data for graphics
   integer :: w_cam  ! Write the fields (0=no; 1=yes)
   integer :: itdt   ! Maximum number of iteractions to update of the time step

   character(200) :: sname !< name of the simulation of the optimizer
   integer :: it     ! counters - iterations
   integer :: optm   ! Coupling with optimization method (0=no; 1=yes)

   integer, parameter :: lid = 100 ! Listing file id
   integer, parameter :: rid = 101 ! Residual file id

   real(8) :: tcpu   ! Total CPU time
   real(8) :: tcpuo  ! Accumulated CPU time before interruption
   real(8) :: tcpu1  ! Starts CPU time
   real(8) :: tcpu2  ! Finishes CPU time

   real(8) :: maxcu  ! Coefficient of convergence for u
   real(8) :: maxcv  ! Coefficient of convergence for v
   real(8) :: maxct  ! Coefficient of convergence for T
   real(8) :: maxcp  ! Coefficient of convergence for p
   real(8) :: maxc   ! Max val among maxcu, maxcv, maxct, maxcp
   real(8) :: maxcc  ! Critical value of maxc

   real(8) :: beta   ! Constant of the UDS/CDS mixing scheme
   real(8) :: dt     ! Time step
   real(8) :: a1     ! Initial step for the geometric progression grid
   real(8) :: rt     ! Radius in the throat of the nozzle [m]

   real(8) :: gamma  ! gamma = Cpo/Cvo in the chamber
   real(8) :: Rg     ! Perfect gas constant
   real(8) :: po     ! Stagnation pressure in the chamber
   real(8) :: T0     ! Stagnation temperature in the chamber
   real(8) :: pr     ! atmospheric pressure at the sea level (Pa)
   real(8) :: go     ! gravitational acceleration at the sea level
   real(8) :: Sg     ! Throttle area
   real(8) :: rcg    ! Throttle radius of curvature

   real(8) :: tolu   ! Tolerance for the MSI method for u, v and T
   real(8) :: tolp   ! Tolerance for the MSI method for p

   real(8) :: fmi    ! Mass flow rate at entrance
   real(8) :: fme    ! Mass flow rate at exit
   real(8) :: Fd     ! Thrust

   real(8) :: fm1D   ! Mass flow rate (Isentropic 1D flow)
   real(8) :: Fd1D   ! Dynamic thrust (Isentropic 1D flow)
   real(8) :: Fpv1D  ! Pressure thrust into vacuum (Isentropic 1D flow)

   real(8) :: norm   ! Norm of the linear systems
   real(8) :: norm1  ! Norm of the linear systems of the first iteraction

   real*8 :: tolerance ! Parada com base na tolerância
   real*8 :: fda_a

   character(len = 3) :: nmbr             ! process number in parallel processing
   character(len = 200) :: folder_input   !< input folder name
   character(len = 200) :: folder_output  !< output folder name
   character(len = 200) :: sim_id         ! Simulation identification
   character(len = 200) :: input_file_parameters ! Input parameters data file
   character(len = 200) :: input_file_constants  ! Input parameters (constants) data file
   character(len = 20) :: day  ! System date
   character(len = 20) :: hour ! System time

   real(8), allocatable, dimension(:) :: x ! X coordinate of the corner ne of each volume (real + ficititious)
   real(8), allocatable, dimension(:) :: y ! Y coordinate of the corner ne of each volume (real + ficititious)
   real(8), allocatable, dimension(:) :: xp ! Centroid X coordinate of each volume (real + ficititious)
   real(8), allocatable, dimension(:) :: yp ! Centroid Y coordinate of each volume (real + ficititious)

   real(8), allocatable, dimension(:) :: xe     ! x_eta at face east of volume P
   real(8), allocatable, dimension(:) :: ye     ! y_eta at face east of volume P
   real(8), allocatable, dimension(:) :: xen    ! x_eta at face north of volume P
   real(8), allocatable, dimension(:) :: yen    ! y_eta at face north of volume P
   real(8), allocatable, dimension(:) :: xk     ! x_csi at face north of volume P
   real(8), allocatable, dimension(:) :: yk     ! y_csi at face north of volume P
   real(8), allocatable, dimension(:) :: xke    ! x_csi at face east of volume P
   real(8), allocatable, dimension(:) :: yke    ! y_csi at face east of volume P
   real(8), allocatable, dimension(:) :: Jp     ! Jacobian at the center of volume P
   real(8), allocatable, dimension(:) :: Je     ! Jacobian at the center of east face of volume P
   real(8), allocatable, dimension(:) :: Jn     ! Jacobian at the center of north face of volume P
   real(8), allocatable, dimension(:) :: Alphae ! (metric) Alpha at the center of east face of volume P
   real(8), allocatable, dimension(:) :: Gamman ! (metric) Gamma at the center of north face of volume P
   real(8), allocatable, dimension(:) :: Betae  ! (metric) Beta  at the center of east face of volume P
   real(8), allocatable, dimension(:) :: Betan  ! (metric) Beta  at the center of north face of volume P

   real(8), allocatable, dimension(:) :: radius ! Radius of northest corner of volume P
   real(8), allocatable, dimension(:) :: re     ! Radius of the center of east face of volume P
   real(8), allocatable, dimension(:) :: rn     ! Radius of the center of north face of volume P
   real(8), allocatable, dimension(:) :: rp     ! Radius of the center of volume P

   real(8), allocatable, dimension(:) :: p      ! Pressure at center of volumes
   real(8), allocatable, dimension(:) :: pa     ! Pressure at center of volumes at a time step before
   real(8), allocatable, dimension(:) :: T      ! Temperature at center of volumes
   real(8), allocatable, dimension(:) :: ro     ! Specific mass (absolute density) at center of volumes
   real(8), allocatable, dimension(:) :: Ta     ! Temperature at the time step before
   real(8), allocatable, dimension(:) :: roa    ! Absolute density at a time step before at center of vol. P
   real(8), allocatable, dimension(:) :: roe    ! Absolute density at east face
   real(8), allocatable, dimension(:) :: ron    ! Absolute density at north face
   real(8), allocatable, dimension(:) :: g      ! g = ro / p

   real(8), allocatable, dimension(:) :: u      ! Cartesian velocity u
   real(8), allocatable, dimension(:) :: v      ! Cartesian velocity v
   real(8), allocatable, dimension(:) :: ua     ! Cartesian velocity of a time step before
   real(8), allocatable, dimension(:) :: va     ! Cartesian velocity of a time step before
   real(8), allocatable, dimension(:) :: ue     ! Cartesian velocity u at center of east face
   real(8), allocatable, dimension(:) :: ve     ! Cartesian velocity v at center of east face
   real(8), allocatable, dimension(:) :: un     ! Cartesian velocity u at center of north face
   real(8), allocatable, dimension(:) :: vn     ! Cartesian velocity v at center of north face
   real(8), allocatable, dimension(:) :: uea    ! Cartesian velocity u at center of east face in the last time step
   real(8), allocatable, dimension(:) :: vea    ! Cartesian velocity v at center of east face in the last time step
   real(8), allocatable, dimension(:) :: una    ! Cartesian velocity u at center of north face in the last time step
   real(8), allocatable, dimension(:) :: vna    ! Cartesian velocity v at center of north face in the last time step

   real(8), allocatable, dimension(:) :: Uce    ! Contravariant velocity U at east face
   real(8), allocatable, dimension(:) :: Vcn    ! Contravariant velocity V at north face

   real(8), allocatable, dimension(:) :: Ucem   ! Contravariant velocity U at east face (previous iteraction)
   real(8), allocatable, dimension(:) :: Vcnm   ! Contravariant velocity V at north face (previous iteraction)

   real(8), allocatable, dimension(:) :: ube !< u over the faces of the east boundary (m/s)
   real(8), allocatable, dimension(:) :: ubw !< u over the faces of the west boundary (m/s)
   real(8), allocatable, dimension(:) :: ubn !< u over the faces of the north boundary (m/s)
   real(8), allocatable, dimension(:) :: ubs !< u over the faces of the south boundary (m/s)

   real(8), allocatable, dimension(:) :: vbe !< v over the faces of the east boundary (m/s)
   real(8), allocatable, dimension(:) :: vbw !< v over the faces of the west boundary (m/s)
   real(8), allocatable, dimension(:) :: vbn !< v over the faces of the north boundary (m/s)
   real(8), allocatable, dimension(:) :: vbs !< v over the faces of the south boundary (m/s)

   real(8), allocatable, dimension(:) :: Tbe !< T over the faces of the east boundary (K)
   real(8), allocatable, dimension(:) :: Tbw !< T over the faces of the west boundary (K)
   real(8), allocatable, dimension(:) :: Tbn !< T over the faces of the north boundary (K)
   real(8), allocatable, dimension(:) :: Tbs !< T over the faces of the south boundary (K)

   real(8), allocatable, dimension(:) :: pbe !< p over the faces of the east boundary (Pa)
   real(8), allocatable, dimension(:) :: pbw !< p over the faces of the west boundary (Pa)
   real(8), allocatable, dimension(:) :: pbn !< p over the faces of the north boundary (Pa)
   real(8), allocatable, dimension(:) :: pbs !< p over the faces of the south boundary (Pa)

   real(8), allocatable, dimension(:) :: Ucbe !< Uc over the faces of the east  boundary (m2/s)
   real(8), allocatable, dimension(:) :: Vcbe !< Vc over the faces of the east  boundary (m2/s)

   real(8), allocatable, dimension(:) :: M      ! Mach number, used for frozen, local equilibrium and non-equilibrium flows

   real(8), allocatable, dimension(:) :: vle    ! Laminar viscosity at center of face east
   real(8), allocatable, dimension(:) :: vln    ! Laminar viscosity at center of face north
   real(8), allocatable, dimension(:) :: vlp    ! Laminar viscosity at center of volume P
   real(8), allocatable, dimension(:) :: cp     ! Specific heat at constant pressure at center o volume P
   real(8), allocatable, dimension(:) :: gcp    ! gcp = gamma = Cp/Cv at center of CV P
   real(8), allocatable, dimension(:) :: kp     ! Thermal conductivity at center of volume P
   real(8), allocatable, dimension(:) :: ke     ! Thermal conductivity at center of face east
   real(8), allocatable, dimension(:) :: kn     ! Thermal conductivity at center of face north

   real(8), allocatable, dimension(:) :: cup    ! Term of deferred correction for u
   real(8), allocatable, dimension(:) :: cvp    ! Term of deferred correction for v
   real(8), allocatable, dimension(:) :: sup    ! Viscous term for u
   real(8), allocatable, dimension(:) :: svp    ! Viscous term for v

   real(8), allocatable, dimension(:) :: pl     ! Pressure correction
   real(8), allocatable, dimension(:) :: due    ! SIMPLEC coefficient for ue
   real(8), allocatable, dimension(:) :: dve    ! SIMPLEC coefficient for ve
   real(8), allocatable, dimension(:) :: dun    ! SIMPLEC coefficient for un
   real(8), allocatable, dimension(:) :: dvn    ! SIMPLEC coefficient for vn
   real(8), allocatable, dimension(:) :: de     ! SIMPLEC coefficient for Uce
   real(8), allocatable, dimension(:) :: dn     ! SIMPLEC coefficient for Vcn
   real(8), allocatable, dimension(:) :: bu     ! Source of the linear system for u
   real(8), allocatable, dimension(:) :: bv     ! Source of the linear system for v
   real(8), allocatable, dimension(:) :: bt     ! Source of the linear system for T
   real(8), allocatable, dimension(:) :: bp     ! Source of the linear system for pl

   real(8), allocatable, dimension(:,:) :: au   ! Coefficients of the linear system for u
   real(8), allocatable, dimension(:,:) :: av   ! Coefficients of the linear system for v
   real(8), allocatable, dimension(:,:) :: at   ! Coefficients of the linear system for T
   real(8), allocatable, dimension(:,:) :: ap   ! Coefficients of the linear system for pl
   real(8), allocatable, dimension(:,:) :: dl9  ! Lower matrix of the MSI method for 9 diagonals
   real(8), allocatable, dimension(:,:) :: du9  ! Upper matrix of the MSI method for 9 diagonals
   real(8), allocatable, dimension(:,:) :: dl5  ! Lower matrix of the MSI method for 5 diagonals
   real(8), allocatable, dimension(:,:) :: du5  ! Upper matrix of the MSI method for 5 diagonals

   real(8), allocatable, dimension(:) :: Twall  ! Wall temperature
   real(8), allocatable, dimension(:) :: Tin    ! Temperature in the entrance
   real(8), allocatable, dimension(:) :: plin   ! Pressure correction in the entrance
   real(8), allocatable, dimension(:) :: uin    ! Cartesian velocity in the entrance
   real(8), allocatable, dimension(:) :: vin    ! Cartesian velocity in the entrance
   real(8), allocatable, dimension(:) :: Mw     ! Mach number in the entrance
   real(8), allocatable, dimension(:) :: pin    ! Pressure in the entrance
   real(8), allocatable, dimension(:) :: pina   ! Pressure in the entrance at a time step before

   real(8), allocatable, dimension(:) :: M1D    ! Mach number distribution. Isentropic flow 1D.
   real(8), allocatable, dimension(:) :: p1D    ! Pressure distribution. Isentropic flow 1D.
   real(8), allocatable, dimension(:) :: T1D    ! Temperature distribution. Isentropic flow 1D.
   real(8), allocatable, dimension(:) :: u1D    ! Velocity distribution. Isentropic flow 1D.

   ! Variables used in the eddy viscosity calculation
   integer :: modtur ! turbulence model option (0=laminar, 1=Baldwin-Lomax)
   real(8), allocatable, dimension(:) :: vtp    ! Eddy viscosity at center of volume P
   real(8), allocatable, dimension(:) :: rona   ! Ro of the previous iteraction level
   real(8), allocatable, dimension(:) :: ro_av  ! Average rho of two iteraction levels
   real(8), allocatable, dimension(:) :: ron_av ! Average rho_n of two iteraction levels
   real(8), allocatable, dimension(:) :: u_av   ! Average u of two iteraction levels
   real(8), allocatable, dimension(:) :: v_av   ! Average v of two iteraction levels

   ! Modifications for chemical reaction flows - new variables
   integer :: pmod ! physical model: 1 - one-species; 2 - one-species (variable properties); 3 - frozen flow; 4 - local equilibrium flow; 5 - non-equilibrium flow
   integer :: cmod ! chemical model (for pmod = 3, 4 or 5)
   real(8) :: OF   ! oxidant/fuel ratio (for pmod = 3, 4 or 5)

   integer :: itimax   ! maximum number of iterations to solve e(j) - dissociation rate for j-reaction
   real(8) :: tol_ej   ! tolerance for e(j) evaluation
   integer :: itemax   ! maximum number of iterations to solve n - the number of mols of products
   real(8) :: tol_n    ! tolerance for n evaluation
   integer :: it_ch    ! number of global iteration to start the reaction evaluation
   integer :: w_ch     ! frequency of global iterations for reaction evaluation
   real(8) :: Tw_cte   ! value for wall temperature (constant)
   integer :: cp_type  ! option for specific heat (1=constant; 2=gases mixture)
   real(8) :: cp_cte   ! value for specific heat (constant)
   integer :: k_type   ! option for themal conductivity (1=constant; 2=gases mixture)
   real(8) :: k_cte    ! value for thermal condutivity (constant)
   integer :: f_type   ! option for friction factor (1=constant; 2=Miller equation (Moody's diagram))
   real(8) :: Darcy    ! value for Darcy's friction factor (constant)
   integer :: v_type   ! option for dynamic viscosity (1=constant; 2=gases mixture)
   real(8) :: visc_cte ! value for dynamic viscosity (constant)
   integer :: h_type   ! option for convection coefficient (1=constant; 2=correlations)
   real(8) :: h_cte    ! value for convection coefficient (constant)
   integer :: rf_type  ! option for recovery factor
   real(8) :: rf_cte   ! value for recovery factor
   real(8) :: eps      ! value of average emissivity: wall-gases
   real(8), parameter :: sigma = 5.67051d-8   ! Stefan-Boltzmann constant

   integer :: Ns       ! Number of species involved for the chosen chemical model

   real(8), allocatable, dimension(:) :: Rgp ! gases mixture constant
   real(8), allocatable, dimension(:,:) :: Yi  ! mass fraction of each chemical species

   character(len = 4), allocatable, dimension(:) :: species ! chemical species involved in the chosen chemical model

   !Modifications for local equilibrium flows - new variables --------------------------------------------------------------------------------------------------

   integer :: Nr       ! Number of chemical reactions for the chosen chemical model
   real(8) :: dpY      ! Pressure variation to estimate gcp
   integer :: app_s    ! Approximation scheme (1=UDS-2/UDS; 2=CDS-2/UDS; 3=Quick/UDS)
   real(8) :: factor   ! Mixture factor for UDS/CDS source term of energy equation

   real(8), allocatable, dimension(:) :: gcpf   ! fgcp = frozen gamma = Cp/Cv at center of CV P (frozen conditions)
   real(8), allocatable, dimension(:) :: entin  ! entropy at the nozzle inlet
   real(8), allocatable, dimension(:) :: ent    ! entropy at each CV P
   real(8), allocatable, dimension(:,:) :: hi     ! enthalpy of each chemical species
   real(8), allocatable, dimension(:,:) :: ej     ! dissociation rate (GIBBS)
   real(8), allocatable, dimension(:,:) :: Yia    ! mass fraction of each chemical species (previous iteration)

   character(len = *), parameter, private :: text_editor = "gedit "

contains

   !****************************************************************************

   ! Subroutine 01

   subroutine get_parameters

      implicit none

      call date_time(day, hour)

      open(10, file = trim(folder_input) // "mach2d-input.dat")

      ! Reading file name from which parameters will be read
      read(10,*) input_file_parameters
      read(10,*) input_file_constants
      read(10,*) pmod   ! Physical model: 1 - one-species; 2 - one-species with variable properties; 3 - frozen flow; 4 - local equilibrium flow; 5 - non-equilibrium flow
      close(10)

      select case (pmod) ! Reading of specific values depending on the physical model

         case(1) ! One-species, constant properties

            !dos = system(text_editor // trim(folder_input) // &
               !trim(adjustl(input_file_parameters)))

            open(10, file = trim(folder_input) // input_file_parameters)
            read(10,*) sim_id ! Simulation identification  (up to 100 characters)
            read(10,*) nx     ! Number of real volumes in the csi direction
            read(10,*) ny     ! Number of real volumes in the eta direction
            read(10,*) kg     ! Kind of grid (1=uniform, 2=Geometric Progression)
            read(10,*) ctd    ! Kind of grid in axial direction (0=uniform; 1=concentred)
            read(10,*) a1     ! Initial step for the geometric progression grid
            read(10,*) coord  ! Coordinate system (1=cylindrical; else=cartesian)
            read(10,*) Rg     ! Perfect gas constant
            read(10,*) gamma  ! gamma = Cpo/Cvo in the chamber (Specific heat ratio)
            read(10,*) po     ! Stagnation pressure in the chamber
            read(10,*) T0     ! Stagnation temperature in the chamber
            read(10,*) pr     ! atmospheric pressure at the sea level (Pa)
            read(10,*) go     ! gravitational acceleration at the sea level
            read(10,*) beta   ! Constant of the UDS/CDS mixing scheme
            read(10,*) modvis ! Viscosity model (0=Euler; 1=Navier-Stokes)
            read(10,*) modtur ! turbulence model option (0=laminar; 1=Baldwin-Lomax)
            read(10,*) dt     ! Initial time step
            read(10,*) itdt   ! Maximum number of iteractions to update of the time step
            read(10,*) maxcc  ! Critical value of maxc
            read(10,*) itmax  ! Maximum number of iteractions of the time evolution
            read(10,*) tolerance ! Parada com base na tolerância
            read(10,*) ittol  ! Number of iterations below the tolerance for stopping criterion
            read(10,*) imax   ! Maximum number of iteractions for the pressure correction
            read(10,*) nitm_u ! Maximum number of iteractions for the MSI method for u, v and T
            read(10,*) tolu   ! Tolerance for the MSI method for u, v and T
            read(10,*) nitm_p ! Maximum number of iteractions for the MSI method for p
            read(10,*) tolp   ! Tolerance for the MSI method for p
            read(10,*) reload ! Upload backup data and continue computation (0=no, 1=yes)
            read(10,*) wbkp   ! Frequency of saving backup data
            read(10,*) wlf    ! Frequency of printing in the listing file
            read(10,*) sem_a  ! Open result files (0=yes; 1=no)
            read(10,*) sem_g  ! Visualize the plot (0=yes; 1=no)
            read(10,*) w_g    ! Frequency of writing data for graphics
            read(10,*) w_cam  ! Write the fields (1=yes; 0=no)
            read(10,*) optm   ! Coupling with optimization method (0=no; 1=yes)
            close(10)

            open(10, file = trim(folder_input) // input_file_constants)
            read(10,*) ccTw    ! ccTw = 0 -> adiabatic;  ccTw = 1 -> prescribed temperature
            read(10,*) Tw_cte  ! prescribed temperature of wall
            read(10,*) cp_cte  ! prescribed constant specific heat
            read(10,*) k_cte   ! prescribed constant thermal conductivity

            select case (modvis)
               case(1) ! Navier-Stokes
                  read(10,*) visc_cte  ! prescribed constant dynamic viscosity
            end select

            close(10)

         case(2) ! One-species, variable properties

            !        dos = system(text_editor // trim(folder_input) // " " // trim(adjustl(input_file_parameters)))

            open(10, file = trim(folder_input) // input_file_parameters)
            read(10,*) sim_id ! Simulation identification  (up to 100 characters)
            read(10,*) nx     ! Number of real volumes in the csi direction
            read(10,*) ny     ! Number of real volumes in the eta direction
            read(10,*) kg     ! Kind of grid (1=uniform; 2=Geometric Progression)
            read(10,*) ctd    ! Kind of grid in axial direction (0=uniform; 1=concentred)
            read(10,*) a1     ! Initial step for the geometric progression grid
            read(10,*) coord  ! Coordinate system (1=cylindrical; else=cartesian)
            read(10,*) specie ! chemical specie (1=H20; 2=AIR)
            read(10,*) Rg     ! Perfect gas constant
            read(10,*) gamma  ! gamma = Cpo/Cvo in the chamber (Specific heat ratio)
            read(10,*) po     ! Stagnation pressure in the chamber
            read(10,*) T0     ! Stagnation temperature in the chamber
            read(10,*) pr     ! atmospheric pressure at the sea level (Pa)
            read(10,*) go     ! gravitational acceleration at the sea level
            read(10,*) beta   ! Constant of the UDS/CDS mixing scheme
            read(10,*) modvis ! Viscosity model (0=Euler; 1=Navier-Stokes)
            read(10,*) modtur ! turbulence model option (0=laminar; 1=Baldwin-Lomax)
            read(10,*) dt     ! Initial time step
            read(10,*) itdt   ! Maximum number of iteractions to update of the time step
            read(10,*) maxcc  ! Critical value of maxc
            read(10,*) itmax  ! Maximum number of iteractions of the time evolution
            read(10,*) tolerance ! Parada com base na tolerância
            read(10,*) ittol  ! Number of iterations below the tolerance for stopping criterion
            read(10,*) imax   ! Maximum number of iteractions for the pressure correction
            read(10,*) nitm_u ! Maximum number of iteractions for the MSI method for u, v and T
            read(10,*) tolu   ! Tolerance for the MSI method for u, v and T
            read(10,*) nitm_p ! Maximum number of iteractions for the MSI method for p
            read(10,*) tolp   ! Tolerance for the MSI method for p
            read(10,*) reload ! Upload backup data and continue computation (0=no; 1=yes)
            read(10,*) wbkp   ! Frequency of saving backup data
            read(10,*) wlf    ! Frequency of printing in the listing file
            read(10,*) sem_a  ! Open result files (0=yes; 1=no)
            read(10,*) sem_g  ! Visualize the plot (0=yes; 1=no)
            read(10,*) w_g    ! Frequency of writing data for graphics
            read(10,*) w_cam  ! Write the fields (0=no; 1=yes)
            read(10,*) optm   ! Coupling with optimization method (0=no; 1=yes)
            close(10)

            ! dos = system(text_editor // trim(folder_input) // " " // trim(adjustl(input_file_constants)))

            open(10, file = trim(folder_input) // input_file_constants)
            read(10,*) ccTw    ! ccTw = 0 -> adiabatic;  ccTw = 1 -> prescribed temperature
            read(10,*) Tw_cte  ! prescribed temperature of wall
            read(10,*) cp_cte  ! prescribed specific heat (initial guess)
            read(10,*) k_cte   ! prescribed thermal conductivity (initial guess)

            select case (modvis)
               case(1) ! Navier-Stokes
                  read(10,*) visc_cte  ! prescribed dynamic viscosity (initial guess)
            end select

            close(10)

         case(3) ! Frozen flows

            ! dos = system(text_editor // trim(folder_input) // trim(adjustl(input_file_parameters)))

            open(10, file = trim(folder_input) // input_file_parameters)
            read(10,*) sim_id ! Simulation identification  (up to 100 characters)
            read(10,*) cmod   ! Chemical model
            read(10,*) nx     ! Number of real volumes in the csi direction
            read(10,*) ny     ! Number of real volumes in the eta direction
            read(10,*) kg     ! Kind of grid (1=uniform; 2=Geometric Progression)
            read(10,*) ctd    ! Kind of grid in axial direction (0=uniform; 1=concentred)
            read(10,*) a1     ! Initial step for the geometric progression grid
            read(10,*) coord  ! Coordinate system (1=cylindrical; else=cartesian)
            read(10,*) Rg     ! Perfect gas constant
            read(10,*) gamma  ! gamma = Cpo/Cvo in the chamber (Specific heat ratio)
            read(10,*) po     ! Stagnation pressure in the chamber
            read(10,*) T0     ! Stagnation temperature in the chamber
            read(10,*) pr     ! atmospheric pressure at the sea level (Pa)
            read(10,*) go     ! gravitational acceleration at the sea level
            read(10,*) beta   ! Constant of the UDS/CDS mixing scheme
            read(10,*) modvis ! Viscosity model (0=Euler; 1=Navier-Stokes)
            read(10,*) modtur ! turbulence model option (0=laminar; 1=Baldwin-Lomax)
            read(10,*) dt     ! Initial time step
            read(10,*) itdt   ! Maximum number of iteractions to update of the time step
            read(10,*) maxcc  ! Critical value of maxc
            read(10,*) itmax  ! Maximum number of iteractions of the time evolution
            read(10,*) tolerance ! Parada com base na tolerância
            read(10,*) ittol  ! Number of iterations below the tolerance for stopping criterion
            read(10,*) imax   ! Maximum number of iteractions for the pressure correction
            read(10,*) nitm_u ! Maximum number of iteractions for the MSI method for u, v and T
            read(10,*) tolu   ! Tolerance for the MSI method for u, v and T
            read(10,*) nitm_p ! Maximum number of iteractions for the MSI method for p
            read(10,*) tolp   ! Tolerance for the MSI method for p
            read(10,*) itimax ! Maximum number of iterations for each chemical reaction
            read(10,*) tol_ej ! Tolerance for each chemical reaction
            read(10,*) itemax ! Maximum number of iterations for the system of chemical reactions
            read(10,*) tol_n  ! Tolerance for the system of chemical reactions
            read(10,*) it_ch  ! Number of global iteration to start the chemical evaluation
            read(10,*) w_ch   ! Frequency of global iterations to chemical evaluation
            read(10,*) OF     ! oxidant/fuel ratio
            read(10,*) reload ! Upload backup data and continue computation (0=no, 1=yes)
            read(10,*) wbkp   ! Frequency of saving backup data
            read(10,*) wlf    ! Frequency of printing in the listing file
            read(10,*) sem_a  ! Open result files (0=yes; 1=no)
            read(10,*) sem_g  ! Visualize the plot (0=yes; 1=no)
            read(10,*) w_g    ! Frequency of writing data for graphics
            read(10,*) w_cam  ! Write the fields (0=no; 1=yes)
            read(10,*) optm   ! Coupling with optimization method (0=no; 1=yes)
            close(10)

            ! dos = system(text_editor // trim(folder_input) // trim(adjustl(input_file_constants)))

            open(10, file = trim(folder_input) // input_file_constants)
            read(10,*) ccTw    ! ccTw = 0 -> adiabatic;  ccTw = 1 -> prescribed temperature
            read(10,*) Tw_cte  ! prescribed temperature of wall
            read(10,*) cp_type ! cp_type = 1 -> constant; cp_type = 2 -> frozen
            read(10,*) cp_cte  ! prescribed constant specific heat
            read(10,*) k_type  ! k_type = 1 -> constant; k_type = 2 -> frozen
            read(10,*) k_cte   ! prescribed constant thermal conductivity

            select case (modvis)
               case(1) ! Navier-Stokes
                  read(10,*) v_type    ! v_type = 1 -> constant; v_type = 2 -> frozen
                  read(10,*) visc_cte  ! prescribed constant dynamic viscosity
            end select

            if (pmod == 4) read(10,*) dpY   ! Pressure variation used to estimate gcp

            close(10)

         case(4) ! Local equilibrium flows

            !        dos = system(text_editor // trim(folder_input) // " " // trim(adjustl(input_file_parameters)))

            open(10, file = trim(folder_input) // input_file_parameters)
            read(10,*) sim_id ! Simulation identification  (up to 100 characters)
            read(10,*) cmod   ! Chemical model
            read(10,*) nx     ! Number of real volumes in the csi direction
            read(10,*) ny     ! Number of real volumes in the eta direction
            read(10,*) kg     ! Kind of grid (1=uniform; 2=Geometric Progression)
            read(10,*) ctd    ! Kind of grid in axial direction (0=uniform; 1=concentred)
            read(10,*) a1     ! Initial step for the geometric progression grid
            read(10,*) coord  ! Coordinate system (1=cylindrical; else=cartesian)
            read(10,*) Rg     ! Perfect gas constant
            read(10,*) gamma  ! gamma = Cpo/Cvo in the chamber (Specific heat ratio)
            read(10,*) po     ! Stagnation pressure in the chamber
            read(10,*) T0     ! Stagnation temperature in the chamber
            read(10,*) pr     ! atmospheric pressure at the sea level (Pa)
            read(10,*) go     ! gravitational acceleration at the sea level
            read(10,*) app_s  ! chosen approximation scheme for the enthalpy contribution at energy equation (1=UDS-2/UDS; 2=CDS-2/UDS; 3=Quick/UDS
            read(10,*) factor ! Constant for the enthalpy contribution mixing scheme (0=UDS; 1=higher order approximation)
            read(10,*) beta   ! Constant of the UDS/CDS mixing scheme
            read(10,*) modvis ! Viscosity model (0=Euler; 1=Navier-Stokes)
            read(10,*) modtur ! turbulence model option (0=laminar; 1=Baldwin-Lomax)
            read(10,*) dt     ! Initial time step
            read(10,*) itdt   ! Maximum number of iteractions to update of the time step
            read(10,*) maxcc  ! Critical value of maxc
            read(10,*) itmax  ! Maximum number of iteractions of the time evolution
            read(10,*) tolerance ! Parada com base na tolerância
            read(10,*) ittol  ! Number of iterations below the tolerance for stopping criterion
            read(10,*) imax   ! Maximum number of iteractions for the pressure correction
            read(10,*) nitm_u ! Maximum number of iteractions for the MSI method for u, v and T
            read(10,*) tolu   ! Tolerance for the MSI method for u, v and T
            read(10,*) nitm_p ! Maximum number of iteractions for the MSI method for p
            read(10,*) tolp   ! Tolerance for the MSI method for p
            read(10,*) itimax ! Maximum number of iterations for each chemical reaction
            read(10,*) tol_ej ! Tolerance for each chemical reaction
            read(10,*) itemax ! Maximum number of iterations for the system of chemical reactions
            read(10,*) tol_n  ! Tolerance for the system of chemical reactions
            read(10,*) it_ch  ! Number of global iteration to start the chemical evaluation
            read(10,*) w_ch   ! Frequency of global iterations to chemical evaluation
            read(10,*) OF     ! oxidant/fuel ratio
            read(10,*) reload ! Upload backup data and continue computation (0=no; 1=yes)
            read(10,*) wbkp   ! Frequency of saving backup data
            read(10,*) wlf    ! Frequency of printing in the listing file
            read(10,*) sem_a  ! open result files (1=no; 0=yes)
            read(10,*) sem_g  ! Vvisualize the plot (0=yes; 1=no)
            read(10,*) w_g    ! Frequency of writing data for graphics
            read(10,*) w_cam  ! Write the fields (0=no; 1=yes)
            read(10,*) optm   ! Coupling with optimization method (0=no; 1=yes)
            close(10)

            open(10, file = trim(folder_input) // input_file_constants)
            read(10,*) ccTw    ! Boundary condition (0=adiabatic; 1=prescribed temperature)
            read(10,*) Tw_cte  ! prescribed temperature of wall
            read(10,*) cp_type ! cp_type = 1 -> constant; cp_type = 2 -> frozen
            read(10,*) cp_cte  ! prescribed constant specific heat
            read(10,*) k_type  ! k_type = 1 -> constant; k_type = 2 -> frozen
            read(10,*) k_cte   ! prescribed constant thermal conductivity

            select case (modvis)
               case(1) ! Navier-Stokes
                  read(10,*) v_type    ! v_type = 1 -> constant; v_type = 2 -> frozen
                  read(10,*) visc_cte  ! prescribed constant dynamic viscosity
            end select

            if (pmod == 4) read(10,*) dpY   ! Pressure variation used to estimate gcp

            close(10)

      end select

      sim_id = adjustl(sim_id)

      ! Number of control volumes in each direction/Total number of control volumes
      nx = nx + 2 ! Real + Fictitious
      ny = ny + 2 ! Real + Fictitious
      nxy = nx*ny

   end subroutine get_parameters

   !****************************************************************************

   ! Subroutine 02

   subroutine write_parameters(fid)
      implicit none
      integer, intent(in) :: fid

      write(fid,*)
      write(fid,*) "Date: ", day
      write(fid,*) "Time: ", hour
      write(fid,*)
      write(fid,*) " ===============================  PARAMETERS  ", &
         "==============================="
      write(fid,*)
      write(fid, "(A20, A)") trim(adjustl(sim_id)), &
         " = sim_id: Simulation identification  (up to 100 characters)"
      write(fid,"(i15, a, A, A)") pmod, &
         " = pmod:   Physical model (1=one-species, constant properties; ", &
         "2=one-species, variable properties; 3=frozen flow; ", &
         "4=local equilibrium flow; 5=non-equilibrium flow)"
      if(pmod>= 3) then
         write(fid, "(i15, a)") cmod, " = cmod:   Chemical model)"
      end if
      write(fid,"(i15, a)") nx-2, &
         " = nx:     Number of real volumes in the csi direction"
      write(fid,"(i15, a)") ny-2, &
         " = ny:     Number of real volumes in the eta direction"
      write(fid,"(i15, a, a)") kg, " = kg:     Kind of grid in ", &
         "radial direction (1=uniform; 2=PG; 3=power function)"
      write(fid,"(i15, a)") ctd, &
         " = ctd:    Kind of grid in axial direction (0=uniform; 1=concentrated)"
      write(fid,"(1pe15.8, a)") a1, &
         " = a1:     Initial step for the geometric progression grid (m)"
      write(fid,"(i15, a)") coord, &
         " = coord:  Coordinate system (1=cylindrical; else=cartesian)"
      write(fid,"(1pe15.8, a)") Rg, &
         " = Rg:     Perfect gas constant (J/(kg.K))"
      write(fid,"(1pe15.8, a)") gamma, &
         " = gamma:  Cpo/Cvo in the chamber (Specific heat ratio)"
      write(fid,"(1pe15.8, a)") po, &
         " = po:     Stagnation pressure in the chamber (Pa)"
      write(fid,"(1pe15.8, a)") T0, &
         " = T0:     Stagnation temperature in the chamber (K)"
      write(fid,"(1pe15.8, a)") pr, &
         " = pr:     Atmospheric pressure at the sea level (Pa)"
      write(fid,"(1pe15.8, a)") go, &
         " = go:     Gravitational acceleration at the sea level (m/s2)"
      if(pmod == 4) then
         write(fid,"(i15, a, a, a)") app_s, " = app_s:  chosen ", &
            "approximation scheme for the enthalpy contribution at energy ", &
            "equation (1=UDS-2/UDS; 2=CDS-2/UDS; 3=Quick/UDS)"
         write(fid,"(1pe15.8, a, a)") factor, &
            " = factor: Constant of the mixing scheme (enthalpy ", &
            "contribution at energy equation) (0=UDS; 1=higher order scheme)"
      end if
      write(fid,"(1pe15.8, a)") beta, &
         " = beta:   Constant of the UDS/CDS mixing scheme"
      write(fid,"(i15, a)") modvis, &
         " = modvis: Viscosity model (0=Euler; 1=Navier-Stokes)"
      write(fid,"(i15, a)") modtur, &
         " = modtur: Turbulence model option (0=laminar; 1=Baldwin-Lomax)"
      write(fid, "(i15, a, a)") ccTw , " = cctw:   Boundary condition for ", &
         "temperature (0=adiabatic; 1=Dirichlet)"
      write(fid,"(1pe15.8, a)") dt, " = dt0:    Initial time step (s)"
      write(fid,"(i15, a)") itdt, &
         " = itdt:   Maximum number of iteractions to update of the time step"
      write(fid,"(1pe15.8, a)") maxcc, &
         " = maxcc:  Critical value of maxc"
      write(fid,"(i15, a)") itmax, &
         " = itmax:  Maximum number of iteractions of the time evolution"
      write(fid,"(1pe15.8, a)") tolerance, &
         " = tolerance: Tolerância utilizada como criterio de parada"
      write(fid,"(i15, a, a)") ittol, " = ittol:  ", &
         "Number of iterations below the tolerance for stopping criterion"
      write(fid,"(i15, a)") imax , &
         " = imax:   Maximum number of iteractions for the pressure correction"
      write(fid,"(i15, a, a)") nitm_u, " = nitm_u: Maximum number of ", &
         "iteractions for the solver for u, v and T"
      write(fid,"(1pe15.8, a)") tolu , &
         " = tolu:   Tolerance for the MSI method for u, v and T"
      write(fid,"(i15, a)") nitm_p, &
         " = nitm_p: Maximum number of iteractions for the MSI method for p"
      write(fid,"(1pe15.8, a)") tolp , &
         " = tolp:   Tolerance for the MSI method for p"
      if(pmod>= 3) then
         write(fid,"(i15, a)") itimax, &
            " = itimax: Maximum number of iterations for each chemical reaction"
         write(fid,"(1pe15.8, a, a)") tol_ej, &
            " = tol_ej: Tolerance for the determination of each chemical ", &
            "reaction dissociation rate"
         write(fid,"(i15, a, a)") itemax, &
            " = itemax: Maximum number of iterations for the evaluation of ", &
            "the number of mols for each product"
         write(fid,"(1pe15.8, a, a)") tol_n, " = tol_n:  Tolerance ", &
            "for the evaluation of the number of mols for each product"
      end if
      write(fid,"(i15, a)") reload, &
         " = reload: Upload backup data and continue computation (0=no; 1=yes)"
      write(fid,"(i15, a)") wbkp, &
         " = wbkp:   Frequency of saving backup data"
      write(fid,"(i15, a)") wlf, &
         " = wlf:    Frequency of printing in the listing file"
      write(fid,"(i15, a)") sem_a, &
         " = sem_a:  Open result files (0=yes; 1=no)"
      write(fid,"(i15, a)") sem_g, &
         " = sem_g:  Visualize the plot (0=yes; 1=no)"
      write(fid,"(i15, a)") w_g, &
         " = w_g:    Frequency of writing data for graphics"
      write(fid,"(i15, a)") w_cam, &
         " = w_cam:  Write the fields (0=no; 1=yes)"
      write(fid,*)

      flush(fid)

   end subroutine write_parameters

   !****************************************************************************

   ! Subroutine 03

   subroutine allocate_variables

      implicit none

      allocate(Tin(ny), pin(ny), pina(ny), plin(ny), uin(ny), vin(ny), Mw(ny))

      allocate(Twall(nx), M1D(nx), p1D(nx), T1D(nx), u1D(nx), g(nxy), vtp(nxy))

      allocate(rona(nxy), ro_av(nxy), ron_av(nxy), u_av(nxy), v_av(nxy))

      allocate(x(nxy), y(nxy), xp(nxy), yp(nxy), xe(nxy), ye(nxy) &
         , xen(nxy), yen(nxy), xk(nxy), yk(nxy), xke(nxy), yke(nxy) &
         , Jp(nxy), Je (nxy), Jn(nxy), alphae(nxy), gamman(nxy), betae(nxy) &
         , betan(nxy), radius(nxy), re(nxy), rn(nxy), rp(nxy), p(nxy) &
         , T(nxy), ro(nxy), roe(nxy), ron(nxy), Uce(nxy), Vcn(nxy) &
         , roa(nxy), vle(nxy), vln(nxy), u(nxy), v(nxy), ua(nxy) &
         , va(nxy), cup(nxy), sup(nxy), vlp(nxy), cvp(nxy), svp(nxy) &
         , cp(nxy), ke(nxy), kn(nxy), Ta(nxy), pl(nxy), bv(nxy), bt(nxy) &
         , bp(nxy), due(nxy), dve(nxy), dun(nxy), dvn(nxy), de(nxy) &
         , dn(nxy), ue(nxy), ve(nxy), un(nxy), vn(nxy), uea(nxy) &
         , vea(nxy), una(nxy), vna(nxy), gcp(nxy), kp(nxy), pa(nxy) &
         , bu(nxy), Ucem(nxy), Vcnm(nxy), M(nxy), Rgp(nxy))

      allocate(ube(ny), ubw(ny), ubn(nx), ubs(nx))
      allocate(vbe(ny), vbw(ny), vbn(nx), vbs(nx))
      allocate(Tbe(ny), Tbw(ny), Tbn(nx), Tbs(nx))
      allocate(pbe(ny), pbw(ny), pbn(nx), pbs(nx))

      allocate(Ucbe(ny), Vcbe(ny))

      allocate(au(nxy,9), av(nxy,9), at(nxy,9), ap(nxy,9), dl9(nxy,5))
      allocate(du9(nxy,4), dl5(nxy,4), du5(nxy,3))

   end subroutine allocate_variables

   !****************************************************************************

   ! Subroutine 04

   subroutine initialize_variables
      implicit none

      folder_input = "./input_files/"
      folder_output = "./output_files/"

      call get_parameters

      call allocate_variables

      ! Initializing variables

      ito = 1

      tcpu = 0.d0

      tcpuo = 0.d0

      x = 0.d0
      y = 0.d0
      xp = 0.d0
      yp = 0.d0

      xe = 0.d0
      ye = 0.d0
      xen = 0.d0
      yen = 0.d0
      xk = 0.d0
      yk = 0.d0
      xke = 0.d0
      yke = 0.d0
      Jp = 0.d0
      Je = 0.d0
      Jn = 0.d0
      alphae = 0.d0
      gamman = 0.d0
      betae = 0.d0
      betan = 0.d0

      radius = 0.d0
      re = 0.d0
      rn = 0.d0
      rp = 0.d0

      p = 0.d0
      T = 0.d0
      ro = 0.d0
      roa = 0.d0
      roe = 0.d0
      ron = 0.d0

      Uce = 0.d0
      Vcn = 0.d0

      cp = 0.d0
      gcp = 0.d0
      kp = 0.d0
      ke = 0.d0
      kn = 0.d0
      vle = 0.d0
      vln = 0.d0
      vlp = 0.d0

      u = 0.d0
      v = 0.d0
      ua = 0.d0
      va = 0.d0

      cup = 0.d0
      cvp = 0.d0
      sup = 0.d0
      svp = 0.d0

      pl = 0.d0
      due = 0.d0
      dve = 0.d0
      dun = 0.d0
      dvn = 0.d0
      de = 0.d0
      dn = 0.d0
      ue = 0.d0
      ve = 0.d0
      un = 0.d0
      vn = 0.d0
      uea = 0.d0
      vea = 0.d0
      una = 0.d0
      vna = 0.d0

      Ucem = 0.d0
      Vcnm = 0.d0

      ube = 0.d0
      ubw = 0.d0
      ubn = 0.d0
      ubs = 0.d0

      vbe = 0.d0
      vbw = 0.d0
      vbn = 0.d0
      vbs = 0.d0

      Tbe = 0.d0
      Tbw = 0.d0
      Tbn = 0.d0
      Tbs = 0.d0

      pbe = 0.d0
      pbw = 0.d0
      pbn = 0.d0
      pbs = 0.d0

      Ucbe = 0.d0
      Vcbe = 0.d0

      au = 0.d0
      av = 0.d0
      at = 0.d0
      ap = 0.d0
      bp = 0.d0
      bu = 0.d0
      bv = 0.d0
      bt = 0.d0

      dl9 = 0.d0
      du9 = 0.d0
      dl5 = 0.d0
      du5 = 0.d0

      Twall = 0.d0
      Tin = 0.d0
      pin = 0.d0
      pina = 0.d0
      plin = 0.d0
      uin = 0.d0
      vin = 0.d0
      Mw = 0.d0

      u1D = 0.d0
      M1D = 0.d0
      p1D = 0.d0
      T1D = 0.d0

      M = 0.d0

      ! Variables used in the eddy viscosity calculation
      vtp = 0.d0
      rona = 0.d0
      ro_av = 0.d0
      ron_av = 0.d0
      u_av = 0.d0
      v_av = 0.d0

   end subroutine initialize_variables

   !****************************************************************************

   ! Subroutine 05

   subroutine date_time(date, time)
      implicit none
      character(len = 20), intent(out) :: date
      character(len = 20), intent(out) :: time
      ! Auxiliary variabless
      integer(4) :: var(8) ! date and time
      character*20 :: vardate, vartime, varzone ! date and time
      character*2 :: aux1, aux2
      character*4 :: aux3
      character*50 :: aux

      call date_and_time(vardate, vartime, varzone, var)

      write(aux,*) var(3)
      aux1 = trim(adjustl(aux))
      write(aux,*) var(2)
      aux2 = trim(adjustl(aux))
      write(aux,*) var(1)
      aux3 = trim(adjustl(aux))
      date = '(' // trim(aux1) // '/' // trim(aux2) // '/' // aux3 // ')'

      write(aux,*) var(5)
      aux1 = trim(adjustl(aux))
      write(aux,*) var(6)
      aux2 = trim(adjustl(aux))
      write(aux,*) var(7)
      aux3 = trim(adjustl(aux))
      time = trim(aux1) // ':' // trim(aux2) // ':' // aux3

   end subroutine date_time

   !============================================================================

end module data
